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ABSTRACT 

Using Constrained Local UniversE Simulations (CLUES) of the formation of the Local Group 
in a cosmological context we investigate the recently highlighted problem that the majority 
of the most massive dark subhaloes of the Milky Way are too dense to host any of its bright 
satellites. In particular, we examine the influence of baryonic processes and find that they 
leave a twofold effect on the relation between the peak of the rotation curve and its posi- 
tion (Vmax and i?i„ax)- Satellites with a large baryon fraction experience adiabatic contraction 
thus decreasing i?niax while leaving Knax more or less unchanged. Subhaloes with smaller 
baryon fractions undergo a decrease in V^ax possibly due to outflows of material. Further- 
more, the situation of finding subhaloes in simulations that lie outside the confidence interval 
for possible hosts of the bright MW dwarf spheroidals, appears to be far more prominent in 
cosmologies with a high a^ normalisation and depends on the mass of the host. We conclude 
that the problem cannot be simply solved by including baryonic processes and hence demands 
further investigations. 

Key words: methods: numerical - A^-body simulations - galaxies: formation - haloes - Local 
Group 



first pointed out by'Klypin et a l.|(|1999) and "Moore et al.|(|1999| ), 
and recently reviewed by|Kravtsov|l 2010^ and|Bullock ( 2010). 



1 INTRODUCTION 

The A Cold Dark Matter (ACDM) model, first explored more than 
two decades ago ( [Davis etal.|1985[ >, has been very successful in ex- 
plaining a multitude of observations at cosmological scales, such as 
anisotropies of Cosmic Microwave Background radiation (CMB) 
[Jarosik & et al.| (e.g. [2011 j and galaxy clustering on large scales 
(e.g. |Cole & et al.|2005| l. However, on smaller, galactic scales the 
tests of the ACDM model are complicated by the baryonic physics 
involved in galaxy formation. Therefore, testing the currently ac- 
cepted concordance model at these scales is necessary in order to 
not only understand the nature of dark matter but also the accuracy 
of the model itself. 

The validity of the ACDM model on galactic scales is still 
being questioned due to the discrepancy between the number of 
observed satellites and the number of predicted dark matter sub- 
haloes. High resolution simulations of galactic-size haloes resolve 
a substantial number of substructures within the virial radius, as 
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The most popular interpretation of this so-called "Missing 
Satellite Problem" requires that the smallest dark matter haloes 
are inefficient at forming stars (e.g. |Bullock|2010]|Kravtsov|2010| ). 
Mechanisms such as early reionization of the intergalactic medium 
and supernovae feedback have been invoked to identify the halo 
mass scale where the galaxy formation starts to be inefficient ( |Bul-| 
[lock et al.|[2000{ |Somerville||2002| [Benson et aL]|2002) , partially 
solving the problem. Furthermore, the detection of satellites is most 
certainly biased because of current detection limits ( [Tollerud et al.[ 
[2008[ [Walsh et al.|2009> . 

There is yet another aspect of the satellite population that 
needs to be addressed: the mismatch between the predicted and in- 
ferred distribution of Knax values at the high-Knax end as recently 
highlighted by ^Boylan-Kolchin et al.[ ( |201 1) , where Vmax measures 
the peak of the rotation curve of subhaloes. Using the Aquarius 
simulations fSpringel et al.|2008| and the Via Lactea II simulation 
( [Diem and et al. 2008 1 they found that the majority of the most mas- 
sive subhaloes (i.e. the high-V^Tiax objects) of the Milky Way are too 
dense to host any of its bright satellites. 

There are a number of ways in which this discrepancy may 
be resolved: the subhalo mass function of the Milky Way could be 
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a statistical anomaly with respect to the ACDM expectations dLiul 
|et al. 2010; Guo et al. 2011 1, or the fundamental assumption that 
the luminosities of the satellites are not monotonically related to 
the mass of the subhaloes does not hold true. 

In response to the claims by [Boylan-Kolchin et al.| pOTTb, 
ILovell et al.| ^2011^ explored the possibility that warm rather than 
cold dark matter can provide a better match to the inferred distri- 
bution of satellite circular velocities. With a power spectrum sup- 
pressed at masses below ~ lO^^M© (corresponding to a warmon 
mass of 2 keV), they found that a warm dark matter model natu- 
rally produces haloes that are less concentrated than their cold dark 
matter counterparts. The attempt to explain the evolution of small 
scale structures in the local universe with a AWDM model was al- 
ready presented in |Tikhonov et al.H2009[ >. However, this is only one 
possible solution to the problem. 

Baryonic processes will most certainly also affect the dark 
matter distribution. [Blumenthal et al.| (1986 1 showed that dissipa- 
tive baryons will lead directly to the adiabatic contraction of the 
halo increasing its central density, thus being a critical ingredient to 
determine subhalo properties. However, the possibility that the in- 
fluence of baryons will lead to a flattening of the dark matter central 
density cusp (through dynamical friction of infalling substructures 
composed of dark matter and baryons) has, for instance, been sug- 
gested by El-Zant et al.| ( [200T] l and further studied in |Romano-Diaz| 
|et al.| ( |2008) . Another way in which the haloes' density can be re- 
duced is through sudden mass outflows that can alter substantially 
the central structure, as suggested by [Navarro et al.| ( (l996[ >. In a re- 
cent work of |Parry et al.| ( |20lT| > this last scenario has been tested by 
following the evolution of one simulated satellite, with promising 
results. The same authors though also showed that the inclusion of 
baryons in simulations does not seem to have any correlation with 
the increase or decrease of the dark matter central density. 

In this Letter we directly address the issue of the Vmax prob- 
lem in ACDM simulations by comparing two identical simulations 
with each other: one that is solely based upon dark matter physics 
and another incorporating all the relevant baryonic physics. These 
simulations form part of the CLUES projecrj in which the initial 
conditions are set by imposing constraints derived from observa- 
tional data of the Local Group. The main feature of using con- 
strained simulations is that it provides a numerical environment that 
closely matches our actual neighborhood. 



2 THE SIMULATIONS 

2.1 Constrained Simulations of the Local Group 



We use the same simulations already presented in Libeskind et al. 
( [20T0t , [Libeskind et"^ ( |20TT] >, [Knebe et al.| \1Q\Q) , and [Knebe 



[et al. ( [2011) and refer the reader to these papers for a more exhaus- 
tive discussion and presentation of these constrained simulations of 
the Local Group that form part of the CLUES project. However, we 
briefly summarize their main properties here for clarity. 

We choose to run our simulations using standard ACDM ini- 
tial conditions, that assume a WMAP3 cosmology I jSpergel et al.[ 
[2007| l, i.e. Q.m = 0.24, 0.^ = 0.042, ^k = 0.76. We use a nor- 
malization of (T8 = 0.75 and a slope of the power spectrum of 
n = 0.95. We used the treePM-SPH MPI code GADGET2 ( [Springel[ 
[2005[ > to simulate the evolution of a cosmological box with side 
length of Lbox ~ 64ft~^Mpc. Within this box we identified (in 



a lower-resolution run utilizing 1024'^ particles) the position of a 
model local group that closely resembles the real Local Group (cf. 
Libeskind et al. 20101. This Local Group has then been re-sampled 
with 64 times higher mass resolution in a region of 2/i~^Mpc about 
its centre giving a nominal resolution equivalent to 4096^ particles 
giving a mass resolution of rn,DM ~ 2.1 x lO^/i~^M0 for the dark 
matter and m^^s = 4.42 x IO^/i^^Mq for the gas particles. For 
more details we refer to the reader to [Gottl6ber et al.[ l|2010). 

For this particular study we further use a gas dynamical SPH 
simulation started with the same initial conditions, in which we ad- 
ditionally follow the feedback and star formation rules of Springel[ 
[& Hernquist (2003) : the interstellar medium (ISM) is modeled as 
a two phase medium composed of hot ambient gas and cold gas 
clouds in pressure equilibrium. The thermodynamic properties of 
the gas are computed in the presence of a uniform but evolving 
ultra-violet cosmic background generated from QSOs and AGNs 
and switched on at z = 6 (Haardt & Madau_19 96). Cooling rates 
are calculated from a mixture of a primordial plasma composition. 
No metal dependent cooling is assumed. Cold gas cloud formation 
by thermal instability, star formation, the evaporation of gas clouds, 
and the heating of ambient gas by supernova driven winds are as- 
sumed to all occur simultaneously. Note that the results presented 
through the paper will only refer to the specific SF/feedback model 
of [Springel & Hemquist| ( [2003| l: other formalisms might lead to 
different conclusions, and will be addressed in a companion paper 
(Di Cintio et al., in preparation). 

In addition we also have at our disposal a dark matter only 
CLUES simulation based upon the WMAP5 cosmology ( [Komatsii] 
[et al.|2009) whose details will be presented in a companion paper; 
here it suffices to know that this simulation has the same formal 
resolution as the WMAP3 one, and it has also been re-simulated 
within a sphere of 2/i~^Mpc radius, i.e. the primary difference be- 
tween the two simulations is merely the cosmology. 



2.2 The (Sub-)Halo Finding 

We used the MPI-l-OpenMP hybrid halo finder AHF (AMIGA halo 
finder) to identify haloes and subhaloes in our simulatior ^ AHF 
is the successor of the MHF halo finder by |Gill et al.| ( [2004 1, and a 
detailed description of its mode of operation is given in the code pa- 
per ( [Knollmann & Knebe[2009[ >. Note that AHF automatically (and 
essentially parameter-free) finds haloes, sub-haloes, sub-subhaloes, 
etc. As the two WMAP3 simulations started with the same initial 
conditions (apart from the baryons) we can match individual sub- 
haloes in the DM only simulation with a "sister" subhalo in the 
SPH run (see [Libeskind etlJ^[20I0[ l. In effect, this cross identi- 
fication pairs subclumps at 2 = that originated from the same 
overdensity in the initial conditions. 



3 RESULTS 



In order to most closely match the results presented by [BoylarTl 
[Kolchin et al.|(20TT| l and not to be contaminated by numerical ef- 
fects we limited the subhaloes used throughout the study to those 
within 300kpc from their respective host and more massive than 
Msub > 2 • lO^M0/i"\ We further stack the data for the two most 
massive hosts representing our MW and M31 galaxies. 

In Fig. [T] we show the relation between -Rmax and Vlnax for 
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Figure 1. The relation between the peak of the rotation curve Vmax and 
its position -Rmax for the WMAP3 simulations: the diamonds are DM only 
subhaloes, the cross represent baryonic SPH subhaloes. The two solid lines 
delimite the Icr confidence interval of the observed bright Milky Way dwarf 
spheroidal galaxies, as in |Boylan-Kolchin et al.|{2bl 11 1. The arrows connect 
the DM-SPH sister pairs found following the matching haloes procedure of 
[Libeskind eTaLKZOTO) . 



the WMAP3 simulation alongside the Icr confidence region of the 
known Milky Way satellites, assuming that the mass density pro- 
file of the subhaloes containing the nine observed dwarf spheroidal 
follows a NFW profile ( [Navarro et al.|1996l l: the two solid lines in 
Fig.[T|(and[2l thus limit the area consistent with the observed half- 
hght radii and masses of these dwarfs, based on the work of |Wolf| 
|etal.|pOTO^ . The diamonds represent the subhaloes in the dark mat- 
ter only simulation while the crosses are the satellite galaxies in the 
SPH run. The lines connect the sister haloes, i.e. those objects that 
could be cross-identified in the two simulations. Please note that 
not all subhaloes could be cross-identified and hence only a cer- 
tain number of them are connected by arrows. The results for the 
WMAP5 (dark matter only) data are presented separately in Fig.|2] 
The results of these plots are quite interesting. First of all we 
notice that massive subhaloes (i.e. high-Vinax objects) appear to be 
outside the observational range only in the case of the WMAP5 
cosmology. This is thus in agreement with the findings of |Boylan^ 
[Kolchin et al.| (J20TTl whose Fig. 2 shows both the Aquarius and 
Via Lactea II simulation data combined. The latter run used a 
erg — 0.74 which is close to our WMAP3 value. Note that for 
this simulation the subhaloes are found only marginally outside 
the observational confidence interval. However, note that the ac- 
tual V^ax values for the subhaloes depend on the host mass and 
Knax.host, respectively (cf. jReed et al.|2005[ |Diemand et al.|2007| 
[Springel et al.|[2008) . Therefore, in order to better compare the 
WMAP5 to the WMAP3 simulation, we scaled the subhaloes' max- 
imum velocities, VZ^%b, by the ratio VZ^l.t/VZ^^st ("ot 
presented here) where the respective values are Knax.MW = 131, 

Knax,M31 = 128 for WMAP3, and K^ax.MW = 178, Kiax.MSl = 

194 for WMAP5 (all in km/s). We find that this re-normalization 
leads to a ~ 30% decrease of the V^^^^^^^ values, bringing them 
into agreement with the WMAP3 results. In that respect, the two 
dark matter only simulations are in fact not too different! 

More importantly, we see in Fig. [T] that the inclusion of bary- 



Figure 2. The same as Fig.fTlbut for the WMAP5 (dark matter only) simu- 
lation. 



onic physics does not solve the problem of the massive and highly 
concentrated dark matter subhaloes. On the contrary subhaloes with 
baryons appear to be down-shifted in the -Rmax-'Knax plane with 
respect to their dark matter counterpart, sometimes even entering 
the regime outside the observational constraints only in the SPH 
run. However, we also find that the lower- Vmax objects seem to be 
shifted in the direction anticipated by [Boylan-Kolchin et al.| ( [2C)l l\ , 
i.e. to the upper left of the plot. There appear to be two competing 
effects moving subhaloes in the -Rmax-Knax plane. 

The six SPH (sister) subhaloes that are outside the confiden- 
tial range have a smaller i?max than their DM only companion: the 
addition of baryons causes a contraction of the halo. This effect 
is also visible for three SPH (sister) subhalo inside the observa- 
tional area and is readily explained by the physical phenomenon of 
adiabatic contraction (Blumenthal et al. 1986, Gnedin et al. 200^. 
We confirm that the baryon fraction /;, = i}i,/Qrn of those sub- 
haloes moving downwards is higher than for the subhaloes shifted 
to the upper left. On average, the baryon fraction of the nine (sis- 
ter) SPH subhaloes, whose -Rmax is reduced with respect to their 
DM counterpart, is /(,//b, cosmic ~ 0.314, while the mean /;, of 
the SPH subhaloes inside the Icr area whose -Rmax increases is 
fb/ fb, cosmic ~ 0.006, i.e. Substantially smaller. The subhaloes 
with high ft experience adiabatic contraction and the majority of 
these objects are the ones with the initial highest i?max - Vmax 
pairs. 

To confirm this last point, we used the CONTRA code jGnedin| 
|et al.|2004l l to calculate the adiabtic contraction of a dark matter 
halo in response to condensation of baryons. Using our numerical 
data as input parameters, we found that adiabatic contraction is ac- 
tually efficient only for those subhaloes with sufficiently high fb, 
as expected: the amount of the 7?max reduction computed this way 
perfectly matches the observed shifts in Fig.fTl 

Instead, for the lower Knax sister subhaloes (with substan- 
tially smaller baryon fractions) the baryonic matter has the capa- 
bility to lower the maximum velocity of the rotation curves, while 
increasing -Rmax- This has already been claimed in previous works 
and may be due to different mechanisms. In particular, we like to 
highlight the mass outflow model of Navarro et al.| ( [T996) : imme- 
diate expulsion of a large fraction of baryonic material during star 
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formation could be the cause of the creation of a central dark matter 
core, which will move the peak of the rotation curve to larger radii. 
This model has been successfully tested by Parry et al. ( 201 1^ who 
followed the formation history of a single stellar dominated satel- 
lite, which undergoes the sequence of events predicted by |NavaiTo| 
|et al.|(T996 ^. Another possible explanation to end up with less con- 
centrated density profiles, is through the mechanism described by 
[Mashchenko et al.] p006 1. A random bulk motion of gas, driven by 
stellar feedback, results in a flattening of the central DM cusps, thus 
leading to DM densities smaller than predicted by pure DM cos- 
mological models. But why is it that those objects with low baryon 
fractions are the ones that require the aforementioned mechanisms? 
Is it that the gas expulsion has already occurred, thereby lowering 
the baryon fraction? Possibly the baryon fraction is only low at red- 
shift z = Q because of mass losses during the subhaloes history? 
Lately, Nickerson et al.| ( [20lT| l explored the effect of several baryon 
loss mechanisms on subhaloes in SPH simulation of a Milky Way 
like galaxy, too: they found that for the subhaloes which ended up 
having (or having had) stars but no gas the most efficient mech- 
anism of baryons removal is exactly the stellar feedback ( |Dekel| 
|& Si lk 1986). Finally, we note that the adiabatic contraction (fol- 
lowing G nedin et al.| ( |2004| ) is ineffective for these subhaloes. We 
will address all these issue of the temporal evolution, mass loss and 
baryon influence in greater detail in a companion paper (Di Cintio 
et al., in preparation). 

We close this Section with a detailed look at the rotation 
curves of the sister haloes in Fig. [3] In each plot the two sister 
objects are presented; the solid and dashed lines represent the cir- 
cular velocity of the DM subhaloes in the dark matter only simu- 
lation and of the (sister) SPH subhalo, respectively. The asterisks 
show the Knax-^max pairs used in Fig.fl] We thus observe adia- 
batic contraction at work: the first three objects (which happen to 
have high baryon fraction) in that plot clearly show the centrally 
peaked total matter distribution in the SPH run. The plot further in- 
dicates that our measurements of the rotation curve and its peak are 
not contaminated by numerical artifacts (e.g. mis-identified halo 
centre, etc.). 



4 CONCLUSIONS 

In this letter we explored the possibility that baryonic processes 
may solve the recently presented problem of "the puzzling dark- 
ness of Milky Way subhaloes" ( [Boylan-Kolchin et aL||2011| l. To 
this extent, we used dark matter only as well as full hydrodynam- 
ical simulations of cosmic structure in the context of the CLUES 
project. We used cosmological parameters determined from both 
the WMAP3 and WMAP5 data. 

Our conclusions are twofold and can be summarized as fol- 
lows: 

• We find that when baryonic physics is included, following 
the feedback and star formation prescriptions of Springel & H ertT] 
Iquist^ (2003) , the problem of having too dense massive subhaloes 
is not solved. Instead, gasdynamical simulations pose new ques- 
tions regarding which mechanisms are responsible for the lower- 
ing of i?max in those subhaloes (while Knax remains more or less 
constant). Adiabatic contraction seems to be a reasonable explana- 
tion, as shown using the modified adiabatic contraction model of 
[Gnedin et al.|(2004l l: this process is effective only for some sub- 
haloes, specifically, for those with a high baryon fraction. For the 
SPH subhaloes with lower baryon fractions at redshift 2 = 0, in- 
stead, we observe a general increase of i?max with respect to their 



DM counterpart, thus meaning that other effects are at work, e.g. 
the model proposed by [Navarro et al.| ([1996} in which a rapid ex- 
pulsion of baryonic mass during star formation causes a reduction 
of the halo concentration, as well as naturally explaining the low 
baryon fraction of these objects. 

• While in the WMAP5 DM only case we find dark matter sub- 
haloes outside the confidence area (calculated following the pre- 
scription given in [Boylan-Kolchin et aLJ ( [2011) ) in the WMAP3 
cosmology we only have one massive subhalo outside this obser- 
vational range. Since the Via Lactea II and Aquarius simulations 
presented in [Boylan-Kolchin et al.[pOTT) are similar cases, we con- 
clude that the cosmology certainly has an influence, too: the higher 
(jg of the WMAP5 scenario eventually led to higher host masses 
which - according to our test - are the most likely reason for the 
higher number of excessively centrally concentrated substructures. 
Note that the latest data from WMAP7 favours ag — 0.807, a value 
between the WMAP3 and WMAP5 results: this could mean that the 
problem is worse than in WMAP3, but not as pronounced as in the 
WMAP5 case. 

An issue neither touched upon by us nor other authors is the 
adequacy of using NFW profiles when calculating the confidence 
interval for possible hosts of the bright MW dwarf spheroidals. It 
is obvious that tidal effects will lead to severe modifications of the 
original NFW density profile subhaloes had upon infall into their 
host jKazantzidis et al.|2004^ . They therefore leave an impact upon 
internal and kinematical properties, respectively ^Lokas etal.|2010| 
[201 1| >, which should be taken into account when using observed 
half-light radii R\/2 and their corresponding masses M1/2 to de- 
fine the confidence interval. Further, [Romano-Diaz et aLJ j2008[ > 
showed that adiabatic contraction makes the dark matter profile al- 
most isothermal. However, the relevance is questionable as material 
will primarily be stripped from the outer regions: [Penarrubia et al.[ 
( ,2008) state that dSphs embedded in NFW haloes are very resilient 
to tidal effects until they are nearly destroyed. This is supported 
by [Navarro et al.[J2010[ > who found that the NFW shape holds rea- 
sonably well even for subhaloes. To roughly gauge the problem, 
we fitted our (SPH) subhaloes to a NFW profile, and observed that 
while some of them are well fitted, there are still objects whose 
density profile cannot be approximated by the simple NFW func- 
tional form. Taking all these considerations into account suggests 
that the NFW profile used to calculate the allowed region is likely 
not the best choice. 

The interpretation of the results presented here clearly de- 
mands a closer investigation of the evolutionary tracks of the satel- 
lites, the actual influence of the SF and feedback model as well as 
an improved calculation of the observational confidence level, ver- 
ifying the applicability of the NFW approach. However, we leave 
this to a companion paper (Di Cintio et al., in preparation) and only 
highlight here that simply the inclusion of baryonic physics does 
not solve the problem; it rather poses new challenges to be explored 
and studied in greater detail. 
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Figure 3. Rotation curves of ten sister pairs of massive subhaloes. In eacli panel the velocity profile of a pair of DM and SPH subhaloes is presented. The 
actual values of Vmax-^max are plotted as asterisks. 
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